Visualization of complex seasonal patterns in time series

Rob J Hyndman

Padova: 23 September 2022

Complex seasonality

Complex seasonality

Complex seasonality

Complex seasonality

Complex seasonal topology

Example: hourly data

Complex seasonality

  • Multiple seasonal periods, not necessarily nested
  • Non-integer seasonality
  • Irregular seasonal topography
  • Seasonality that depends on covariates
  • Complex seasonal topology
  • How to effectively visualise the underlying seasonalities?
  • How to decompose such time series into trend and multiple season components?

Visualizing complex seasonalities

Granularity plots

Granularity plots

Granularity plots

Granularity plots

Granularity plots

Granularity plots

Granularity plots

Granularity plots

Faceted granularity plots

Faceted granularity plots

Faceted granularity plots

Faceted granularity plots

Faceted granularity plots

Granularities

Available cyclic granularities for half-hourly data:

 [1] "hhour_hour"         "hhour_day"          "hhour_week"        
 [4] "hhour_fortnight"    "hhour_month"        "hhour_quarter"     
 [7] "hhour_semester"     "hhour_year"         "hour_day"          
[10] "hour_week"          "hour_fortnight"     "hour_month"        
[13] "hour_quarter"       "hour_semester"      "hour_year"         
[16] "day_week"           "day_fortnight"      "day_month"         
[19] "day_quarter"        "day_semester"       "day_year"          
[22] "week_fortnight"     "week_month"         "week_quarter"      
[25] "week_semester"      "week_year"          "fortnight_month"   
[28] "fortnight_quarter"  "fortnight_semester" "fortnight_year"    
[31] "month_quarter"      "month_semester"     "month_year"        
[34] "quarter_semester"   "quarter_year"       "semester_year"

Plot options:

  • raw data or distributional summary on y-axis
  • granularity on x-axis
  • optional granularity as facet

Single granularity plots

  • See distributional differences between categories.
  • Measure the effectiveness of a plot as the average Jenson-Shannon divergence between categories (or between neighbouring categories).
  • Significant differences can be measured using permutation tests or \chi^2 tests.
  • Users can be guided to view the most effective plots.

Faceted granularity plots

  • Omit combinations with empty or near-empty intersections (“clashes”). e.g., day_year \times month_year
  • Prioritise within panel differences over between panel differences
  • Rank combinations using normalized Median Maximum Pairwise Distances.

MSTL

 

  • Kasun Bandara, Rob J Hyndman, Christoph Bergmeir (2022) MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns. International J Operational Research, to appear. robjhyndman.com/publications/mstl/
  • For multiple integer seasonal periods with additive components
  • Implemented in R packages forecast and fable.

MSTL

velec |>
  model(STL(Demand)) |>
  components() |>
  autoplot()

MSTL

MSTL

y_t = T_t + \sum_{i=1}^I S_t^{(i)} + R_t

y_t= observation at time t
T_t= smooth trend component
S_t^{(i)}= seasonal component i
i = 1,\dots,I
R_t= remainder component

Estimation

Components updated iteratively.

MSTL

# X: time series object
# periods: vector of seasonal periods in increasing order
# s.window: seasonal window values
# iterate: number of  STL iterations

seasonality <- matrix(0, nrow = nrow(X), ncol = length(periods))
deseas <- X
for (j in 1:iterate) {
  for (i in 1:length(periods)) {
    deseas <- deseas + seasonality[, i]
    fit <- model(
      STL(deseas ~ season(period = periods[i], window = s.window[i]))
    ) |>
      components()
    seasonality[, i] <- fit$season
    deseas <- deseas - seasonality[, i]
  }
}
trend <- fit$trend
remainder <- deseas - trend
return(trend, remainder, seasonality)

MSTL

fable syntax

tsibble |>
  model(STL(variable) ~ season(period = a, window = b) +
                        season(period = c, window = d))



forecast syntax

vector |>
  msts(seasonal.periods = c(a, c)) |>
  mstl(s.window = c(b, d))

STR

 

  • Alex Dokumentov and Rob J Hyndman (2022) STR: Seasonal-Trend decomposition using Regression. INFORMS Journal on Data Science, to appear. robjhyndman.com/publications/str/

  • Implemented in R package stR.

STR

y_{t} = T_{t} + \sum_{i=1}^{I} S^{(i)}_{t} + \sum_{p=1}^P \phi_{p,t} z_{t,p} + R_{t}

T_t= smooth trend component
S_t^{(i)}= seasonal component i (possibly complex topology)
z_{p,t}= covariate with coefficient \phi_{p,t} (possibly time-varying)
R_t= remainder component

Estimation

Components estimated using penalized MLE

Smoothness via difference operators

Smooth trend obtained by requiring \Delta_2 T_t \sim \text{NID}(0,\sigma_L^2)

  • \Delta_2 = (1-B)^2 where B= backshift operator
  • \sigma_L controls smoothness

f(\bm{D}_\ell \bm{\ell}) \propto \exp\left\{-\frac{1}{2}\big\|\bm{D}_\ell \bm{\ell} / \sigma_L\big\|_{L_2}^2\right\}

  • \bm{\ell} = \langle T_{t} \rangle_{t=1}^{n}
  • \bm{D}_\ell= 2nd difference operator matrix: \bm{D}_\ell\bm{\ell} = \langle\Delta^2 T_{t}\rangle_{t=3}^n

Smooth 2D seasonal surfaces

  • m_i= number of “seasons” in S^{(i)}_{t}.
  • S^{(i)}_{k,t}= 2d season (k=1,\dots,m_i;t=1,\dots,n)
  • \sum\limits_k S^{(i)}_{k,t} = 0 for each t.

Smooth 2D seasonal surfaces

  • \bm{S}^{(i)} = [S_{k,t}^{(i)}] the ith seasonal surface matrix
  • \bm{s}_i = \text{vec}(\bm{S}_i)= the ith seasonal surface in vector form

Smoothness in time t direction:

\begin{align*} \bm{D}_{tt,i} \bm{s}_i &= \langle \Delta^2_{t} \bm{S}^{(i)}_{k,t} \rangle \sim \text{NID}(\bm{0},\sigma_{i}^2 \bm{\Sigma}_{i})\\ f(\bm{s}_i) &\propto \exp\Big\{-\frac{1}{2}\big\|\ \bm{D}_{tt,i}\bm{s}_i / \sigma_i\big\|_{L_2}^2\Big\} \end{align*}

Analogous difference matrices \bm{D}_{kk,i} and \bm{D}_{kt,i} ensure smoothness in season and time-season directions.

Gaussian remainders

  • R_{t} \sim \text{NID}(0,\sigma_R^2).
  • \bm{y} = [y_1,\dots,y_n]'= vector of observations
  • \bm{Z}=[z_{t,p}]= covariate matrix with coefficient \bm{\Phi} = [\phi_{p,t}]
  • \bm{Q}_i= matrix that extracts \langle S^{(i)}_{\kappa(t),t} \rangle_{t=1}^{n} from \bm{s}_i.
  • Residuals: \bm{r} = \bm{y} - \sum_i\bm{Q}_i\bm{s}_i -\bm{\ell} - \bm{Z}\bm{\Phi} have density f(\bm{r}) \propto \exp\Big\{-\frac{1}{2}\big\|\bm{r}/\sigma_R\big\|_{L_2}^2\Big\},

MLE for STR

Minimize wrt \bm{\Phi}, \bm{\ell} and \bm{s}_i:

\begin{align*} -\log \mathcal{L} &= \frac{1}{2\sigma_R} \Bigg\{ \Big\| \bm{y}- \sum_{i=1}^I \bm{Q}_i\bm{s}_i - \bm{\ell} - \bm{Z}\bm{\Phi} \Big\|_{L_2}^2 + \lambda_\ell\Big\|\bm{D}_\ell \bm{\ell}\Big\|_{L_2}^2 \\ & \hspace*{1cm} + \sum_{i=1}^{I}\left( \left\|\lambda_{tt,i} \bm{D}_{tt,i} \bm{s}_i \right\|_{L_2}^2 + \left\|\lambda_{st,i} \bm{D}_{st,i} \bm{s}_i \right\|_{L_2}^2 + \left\|\lambda_{ss,i} \bm{D}_{ss,i} \bm{s}_i \right\|_{L_2}^2 \right) \Bigg\} \end{align*}

Equivalent to linear model

\bm{y}_{+} = \bm{X}\bm{\beta} + \bm{\varepsilon}

  • \bm{y}_{+} = [\bm{y}',~ \bm{0}']'
  • \bm{\varepsilon} \sim N(\bm{0},\sigma_R^2\bm{I})

\bm{X} = \begin{bmatrix} \bm{Q}_1 & \dots & \bm{Q}_I & \bm{I}_n & \bm{Z} \\ \lambda_{tt,1} \bm{D}_{tt,1} & \dots & 0 & 0 & 0 \\ \lambda_{st,1} \bm{D}_{st,1} & \dots & 0 & 0 & 0 \\ \lambda_{ss,1} \bm{D}_{ss,1} & \dots & 0 & 0 & 0 \\ 0 & \ddots & 0 & 0 & 0 \\ 0 & \dots & \lambda_{tt,I} \bm{D}_{tt,I} & 0 & 0 \\ 0 & \dots & \lambda_{st,I} \bm{D}_{st,I} & 0 & 0 \\ 0 & \dots & \lambda_{ss,I} \bm{D}_{ss,I} & 0 & 0 \\ 0 & \dots & 0 & \lambda_\ell \bm{D}_{tt} & 0 \end{bmatrix}

STR

Three seasonal components, quadratic temperature regressors

STR outliers


For more information

Slides: robjhyndman.com/seminars/padova2022

Find me at:

robjhyndman.com
@robjhyndman
@robjhyndman
rob.hyndman@monash.edu